Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a danymi pomocniczymi stosowanymi w modelowaniu dazymetrycznym

1. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a punktami adresowymi.

maxpop = max(c(out_model_punkty_adresowe$TOT, out_model_punkty_adresowe$EST_POP))

p1 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = TOT)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) + 
  labs(title = "Ryc.1 Mapa oryginalnych wartości liczby ludności") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

p2 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = EST_POP)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop)) + 
  labs(title = "Ryc.2 Mapa estymowanych wartości liczby ludności") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p1)

print(p2)

p3 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = RES)) +
  scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(out_model_punkty_adresowe$RES), max(out_model_punkty_adresowe$RES))) + 
  labs(title = "Ryc.3 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p3)

out_model_punkty_adresowe$RES_reclass = case_when(
    out_model_punkty_adresowe$RES < 0 ~ "bledy ujemne",
    out_model_punkty_adresowe$RES == 0 ~ "brak bledu",
    out_model_punkty_adresowe$RES > 0 ~ "bledy dodatnie",
  )

p4 <- ggplot(data = out_model_punkty_adresowe) +
  geom_sf(aes(fill = RES_reclass)) +
  scale_fill_manual(name = "Typ błędu", 
                    values = c("bledy dodatnie" = "red", "brak bledu" = "white", "bledy ujemne" = "blue")) +
  labs(title = "Ryc.4 Mapa zreklasyfikowanych błędów modelu") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

print(p4)

meanerror = mean(out_model_punkty_adresowe$RES)
rmse = sqrt(mean((out_model_punkty_adresowe$RES)^2))

print(paste("średni błąd estymacji modelu wynosi", round(meanerror,2), "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))

## [1] "średni błąd estymacji modelu wynosi 3.67 natomiast pierwiastek średniego błędu kwadratowego wynosi 30.44"

2. Wykorzystanie analizy regresji do określenia zależności między ogólną liczbą ludności a liczbą mieszkań.

2.1 Liniowy model zależności między liczbą ludności a liczbą mieszkań w siatce 1km

# wykres rozrzutu
ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

model_liniowy <- lm(TOT ~ N_BUDYNKI, data = pop1km)

par(mfrow = c(2, 2))
plot(model_liniowy, which = 1)
plot(model_liniowy, which = 2)
plot(model_liniowy, which = 3)
plot(model_liniowy, which = 4)
mtext("Ryc. 5 Wykresy diagnostyczne modelu liniowego", outer = TRUE, line = -1.5, cex = 1)

W jakim stopniu liczba ludności (TOT) jest wyjaśniana przez liczbę mieszkań?

2.2. Jakie są statystyki reszt? (residuals)

Ryc.6 Interaktywne wykresy reszt

## Statystyki

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -98.6520  -2.5883  -0.7094   0.0000  -0.5883 174.7735

## Średnia:  3.738246e-16

3. Ocena dopasowania modelu (wykresy diagnostyczne, identyfikacja wartości odstających).

Wartości odstające można zidentyfikować wykorzystując wykresy diagnostyczne, a w szczególności wykres Residuals vs. Leverage.

p <- autoplot(model_liniowy) + theme_classic()

gridExtra::grid.arrange(grobs = p@plots, top = "Ryc. 7 Wykresy diagnostyczne modelu liniowego")

#dopasowanie modelu
smr_model_lm <-summary(model_liniowy)

c( R2 = smr_model_lm$r.squared, R2_adj = smr_model_lm$adj.r.squared)

##        R2    R2_adj 
## 0.7019419 0.6990481

~70% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI.

# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy, level = 0.99)

##                 0.5 %   99.5 %
## (Intercept) -7.022728 8.441527
## N_BUDYNKI    3.275627 4.603257

outlierTest(model_liniowy)

##     rstudent unadjusted p-value Bonferroni p
## 21  9.742598         3.0248e-16   3.1760e-14
## 25  6.813338         6.7995e-10   7.1395e-08
## 34 -5.176442         1.1389e-06   1.1958e-04

ggplot(pop1km, aes(x=TOT, y=N_BUDYNKI)) + 
  geom_point() +
  geom_point(data = pop1km[c(21, 25, 34),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
  labs(title = "Ryc.8 Wykres rozrzutu z oznaczonymi punktami o wartościach odstających") +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

meanerror = mean(model_liniowy$residuals)
rmse = sqrt(mean((model_liniowy$residuals)^2))

print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))

## [1] "średni błąd estymacji modelu wynosi 3.73824648485021e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 26.25"

PYTANIA Czego dowiadujemy się z wykresów diagnostycznych? Czy model jest dobrze dopasowany? Czy spełnione są założenia regresji liniowej?

4. Model po usunięciu wartości odstających

4.1 Model regresji

# Usunięcie wierszy 21, 25 i 34
pop1km_n <- pop1km[-c(21, 25, 34), ]

# wykres rozrzutu
ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + geom_point() +
  labs(title = "Ryc.9 Wykres rozrzutu po usunięciu punktów z wartościami odstającymi") +
  coord_fixed(ratio = 1) + 
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# model liniowy
model_liniowy_n <- lm(TOT ~ N_BUDYNKI, data = pop1km_n)

p1 <- autoplot(model_liniowy_n) + theme_classic()

gridExtra::grid.arrange(grobs = p1@plots, top = str_wrap("Ryc. 10 Wykresy diagnostyczne modelu liniowego po usunięciu punktów z wartościami odstającymi"))

4.2. Ocena modelu

#dopasowanie modelu
smr_model_lm_n <-summary(model_liniowy_n)

c( R2 = smr_model_lm_n$r.squared, R2_adj = smr_model_lm_n$adj.r.squared)

##        R2    R2_adj 
## 0.8914342 0.8903486

~89% zmienności zmiennej TOT jest wyjaśniona przez zmienność zmiennej N_BUDYNKI.

# Określenie dowolnego przedziału ufności dla współczynników modelu
confint(model_liniowy_n, level = 0.99)

##                 0.5 %   99.5 %
## (Intercept) -2.413594 3.198192
## N_BUDYNKI    3.198037 3.843294

outlierTest(model_liniowy_n)

##     rstudent unadjusted p-value Bonferroni p
## 65 -6.324628         7.3790e-09   7.5266e-07
## 33  4.742867         7.0854e-06   7.2271e-04

ggplot(pop1km_n, aes(x=TOT, y=N_BUDYNKI)) + 
  geom_point() +
  geom_point(data = pop1km_n[c(65, 33),], aes(x=TOT, y=N_BUDYNKI), color = "red", size = 4) +
  labs(title = str_wrap("Ryc.11 Wykres rozrzutu z oznaczonymi punktami o wartościach odstających (po usunięciu punktów)")) +
  xlab("Liczba ludności") + 
  ylab("Liczba mieszkań") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#WIZUALIZACJA
par(mfrow = c(2, 4))
plot(model_liniowy)
plot(model_liniowy_n)
mtext(str_wrap("Ryc. 12 Wykresy diagnostyczne modelu liniowego (pierwszy wiersz - przed usunięciem, drugi wiersz - po)"), adj = 1, line = 20.5, font = 2, cex = 0.8)

meanerror = mean(model_liniowy_n$residuals)
rmse = sqrt(mean((model_liniowy_n$residuals)^2))

print(paste("średni błąd estymacji modelu wynosi", meanerror, "natomiast pierwiastek średniego błędu kwadratowego wynosi", round(rmse,2)))

## [1] "średni błąd estymacji modelu wynosi 2.30805383656037e-16 natomiast pierwiastek średniego błędu kwadratowego wynosi 9.14"

5. Wizualizacja wyników modelu na mapie

5.1. Mapa pokazująca oryginalne wartości liczby ludności (TOT)

maxpop1 = max(c(pop1km$TOT, pop1km$EST_POP))

ggplot() +
  geom_sf(data = pop1km, aes(fill = TOT)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop1)) +
  labs(fill = "Liczba ludności", title = "Ryc.13 Mapa oryginalnych wartości liczby ludności") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

5.2. Mapa pokazującą estymowaną liczbę ludności (EST_POP)

# Przewidywanie EST_POP na podstawie modelu liniowego
pop1km$EST_POP <- predict(model_liniowy, newdata = pop1km)

# Tworzenie mapy z estymowaną liczbą ludności (EST_POP)
ggplot() +
  geom_sf(data = pop1km, aes(fill = EST_POP)) +
  scale_fill_gradient2(name = "Liczba ludności", low = "white", high = "red", limits = c(0,maxpop1)) +
  labs(fill = "Estymowana liczba ludności", title = "Ryc.14 Mapa estymowanych wartości liczby ludności") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

5.3. Mapa reszt (różnic między wartością obserwowaną TOT oraz EST_POP)

# Obliczenie reszt
pop1km$Reszty <- pop1km$TOT - pop1km$EST_POP

# Tworzenie mapy reszt
ggplot() +
  geom_sf(data = pop1km, aes(fill = Reszty)) +
  scale_fill_gradient2(name = "Błąd estymacji", low = "blue", high = "red", limits = c(min(pop1km$Reszty), max(pop1km$Reszty))) +
  labs(fill = "Reszty", title = "Ryc.15 Mapa różnic między wartością obserwowaną liczby ludności a jej estymacją") +
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

5.4. Mapa reszt przeklasyfikowaną do 3 klas (błędy ujemne, brak błędu, błędy dodatnie)

# Przeklasyfikowanie reszt do trzech klas
pop1km$Klasy_reszt = case_when(
    pop1km$Reszty < 0 ~ "Niedoszacowane",
    pop1km$Reszty == 0 ~ "Brak błędu",
    pop1km$Reszty > 0 ~ "Przeszacowane",
  )

# Tworzenie mapy przeklasyfikowanych reszt
ggplot() +
  geom_sf(data = pop1km, aes(fill = Klasy_reszt)) +
  scale_fill_manual(name = "Typ błędu", values = c("Przeszacowane" = "red", "Brak błędu" = "white", "Niedoszacowane" = "blue")) +
  labs(title = "Ryc.16 Mapa zreklasyfikowanych błędów modelu") + 
  theme_bw() + theme(
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank())

6. Statystyczna i przestrzenna analiza rokładu błędów modelu regresji na podstawie map oraz wyników modelu

## Statystyki reszt:

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -98.6520  -2.5883  -0.7094   0.0000  -0.5883 174.7735

## Średnia reszt:  -5.399296e-16

## Odchylenie standardowe reszt:  26.37129

# Histogram reszt
ggplot(pop1km, aes(x=Reszty)) +
  geom_histogram(bins=30, fill="blue", color="black") +
  labs(title = "Ryc.17 Histogram reszt") +
  xlab("Reszty") +
  ylab("Częstość") +
  theme_minimal()

7. Mapa rozmieszczenia ludności w siatce 100m, dane pomocnicze - liczba mieszkań.

leaflet(data = pop_100m1, width = "100%", height = "600px") %>%
  addTiles() %>%
  addPolygons(
    data = pop_100m1,
    fillColor = ~pal(TOT),
    group = "Liczba ludności",
    color = "#BDBDC3",
    weight = 1,
    opacity = 1,
    fillOpacity = 1,
    highlightOptions = highlightOptions(
      weight = 2,
      color = "#666",
      fillOpacity = 1,
      bringToFront = TRUE
    ),
    label = ~paste("Liczba ludności:", TOT),
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  addLegend(
    pal = pal,
    values = ~TOT,
    opacity = 1,
    title = "Liczba ludności",
    position = "bottomright"
  ) %>%
  addLayersControl(
    baseGroups = c("Podkład mapy"),
    overlayGroups = c("Liczba ludności"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  setView(lng = mean(st_coordinates(pop_100m1)[,1]), 
          lat = mean(st_coordinates(pop_100m1)[,2]), 
          zoom = 12)

8. Porównanie wyników obu modeli

Wykonanie analizy regresji liniowej określającej zależność między ogólną liczbą ludności a liczbą mieszkań daje dokładniejszy model, niż zależność między ogólną liczbą ludności a punktami adresowymi. Pomimo że w drugim modelu nie występują pola pozbawione błędu, to średni błąd wciąż jest znacznie mniejszy niż w przypadku pierwszego modelu. Z wykresów dowiadujemy się, że występują wartości odstające, które w znaczny sposób utrudniają odpowiednie dopasowanie modelu do danych, jak i wpływają na wyniki estymacji. Po usunięciu wartości odstających modele są lepiej dopasowane. W resztach są widoczne wzorce, na Q-Q Plot widoczne są odchylenia, reszty wykazują zmienną wariancję, a także obecne są wartości wpływowe. Te czynniki sprawiają, że nawet po usunięciu wartości odstających założenia regresji liniowej nie są spełnione wcale, bądź są spełnione tylko w części.